Surface superconductivity in multilayered rhombohedral graphene: Supercurrent 
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The supercurrent for the surface superconductivity of a fiat-band multilayered rhombohedral 
graphene is calculated. Despite the absence of dispersion of the excitation spectrum, the supercur- 
rent is finite. The critical current is proportional to the zero-temperature superconducting gap, i.e., 
to the superconducting critical temperature and to the size of the fiat band in the momentum space. 
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I. INTRODUCTION 

Fermionic systems with dispersionless branches of ex- 
citation spectrum (flat bands) have quite unusual proper- 
ties; nowadays they attract lots of research interest. Flat 
bands were predicted in many condensed matter systems, 
see for example Refs. 1-4 . In some cases the flat bands are 
protected by topology in momentum space; they emerge 
on the surfaces of gapless topological matter— such as 
surfaces of nodal superconductors^, graphene edges^, 
surfaces of multilayered graphene structures^—, and in 
the cores of quantized vortices in topological superfluids 
and superconductor o 5 ' 1:L i 12 . 

The singular density of states (DOS) associated with 
the dispersionless spectrum may essentially enhance the 
transition temperature opening a new route to room- 
temperature superconductivity. The corresponding crit- 
ical temperature depends linearly on the pairing inter- 
action strength and can be thus considerably higher 
than the usual exponentially small critical temperature 
in the bulki&H. It was shown ir&H that the flat band 
that appears on the surface of multilayered rhombohe- 
dral graphene is especially favorable for surface supercon- 
ductivity. Formation of surface superconductivity is en- 
hanced already for a system having N > 3 layers, where 
the normal-state spectrum has a slow power-law disper- 
sion £ p cx IpI^ as a function of the in-plane momentum 
p. The DOS i/(f p ) cx £ ( p~ N)/N has a singularity at zero 
energy which results in a drastic enhancement of the crit- 
ical temperature. 

Absence of dispersion in a flat band raises the ques- 
tions of superconducting velocity and of the supercur- 
rent: Can they be nonzero and, if they can, what is then 
the magnitude of the critical current? In this Letter we 
address the problem of supercurrent associated with the 
surface superconductivity in the flat-band multilayered 
rhombohedral graphene. Based on the model employed 
in Refii^ for description of the surface superconductivity 
we calculate the supercurrent as a response to a small 
gradient of the order parameter phase using an approach 
similar to that used for calculations of the supercurrent 
in a single layer of graphene 1 ^. We demonstrate that the 
supercurrent is finite; the critical current is proportional 
to the superconducting zero-temperature gap, i.e., to the 



critical temperature, and to the radius of the flat band 
in the momentum space. Being produced by the surface 
superconductivity, the total current through the sample 
is independent of the sample thickness. 



II. THE MODEL 

As in Ref^ we consider multilayered graphene struc- 
ture of N layers in the discrete representation with 
respect to interlayer coupling. For simplicity we 
choose the rhombohedral stacking configuration consid- 
ered in££ ~ 10 i 13 and assume that the most important are 
hoppings between the atoms belonging to different sub- 
lattices parameterized by a single hopping energy t. More 
general form of the multilayered Hamiltonian can be 
found in Refsj 15 ' . In the superconducting case the 
Hamiltonian has the form of a matrix in the Nambu 
space. The Bogoliubov-de Gennes (BdG) equations are 
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where the sum runs over the layers. The normal-state 
Hamiltonian^ 

Hij = v F (& ■ p)S hJ - ta + 5ij + i - taS^j-i , (1) 

<t = (<j x , & y ), cr± = (a x ±ia y )/2 J and Vi are matrices 
and spinors in the pseudo-spin space associated with two 
sublattices. This Hamiltonian acts on the envelope func- 
tion of the in-plane momentum p taken near one of the 
Dirac points, i.e., for |p| <C h/a where a is the interatomic 
distance within a layer; vf = 3toa/2h where to is the the 
hopping energy between nearest-neighbor atoms belong- 
ing to different sublattices on a layer. The particle-like, 
Ui, and hole- like, wave functions near the Dirac point 
are coupled via the superconducting order parameter A, 
that can appear in the presence of a pairing interaction. 
Here we do not specify the nature of the pairing. It can be 
due to either electron-phonon interaction or other pair- 
ing interactions that have been suggested as a source for 
intrinsic superconductivity in graphene, see Refs.— . The 
excitation energy for particles and holes is measured up- 
wards or downwards, respectively, from the Fermi level 
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which can be shifted with respect to the Dirac point due 
to doping. Here we assume that the shift is the same on 
all layers. The order parameter and the Fermi level shift 
/i are scalars in the pseudo-spin space. We assume that 
A and fi are much smaller than the inter-layer coupling 
energy t > 0, which in turn is t -C to. Usually, t ~ 0.1 to 
where to ~ 3 eV— . 

We decompose the wave function 



This gives 

d [(A + )^A+ + (i-) f i-] = l - v 2 FP 2 /t 2 . (8) 

A finite order parameter A couples the particle and 
hole channels at the outermost layers, i — 1 and i = N, 

faVF(px - ip v )&T - Tafii&f = Ea~l - Aaf , (9) 
f 3 v F (p x + ip y )a^f - f 3 [i N a N = Ea N - Aa N . (10) 
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into the spinor functions localized at each sublattice 
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We introduce matrices and vectors in the Nambu space 
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The BdG equations take the form 

f 3 [vf(Px - iPy)Un ~ t&n-l ~ ^ & n] = E& t> " ¥= 1 , (3) 



The boundary conditions (|9]), (|T0|) select p z and deter- 
mine 2N particle and hole branches of the energy spec- 
trum. Looking for the branches that belong to the sur- 
face states with energies of the order of A and n, we solve 
these equations for E <C t. 



III. SUPERCURRENT 



The operator of current along a layer couples 
the states at different sublattices, vl p (n)&u 7 ^p(n) + 
wJ iP (n)(TUy lP (n). For example, the x component of cur- 
rent at layer n is 



?3 [vf(Px + iPy)®! 



ta 
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Ea n , n^N, (4) j x (n) = -ev F V [&+i n (p)a (p) + a * (p)a+ n (p)] 



where p is the momentum operator. In Eqs. ([3]) and (|4]) 
we assume that A„ 7^ only at the outermost layers, 
while A„ = for n ^ 1,N. The arguments supporting 
this assumption are given in RefJ^; it was shown that 
the order parameter quickly decays as a function of the 
distance from the surface. We also neglect A„ as com- 
pared to t in Eqs. © and (01 for n — N and n = 1, 
respectively, as they lead to higher-order corrections in 
A/t. The particle and hole channels are thus decou- 
pled if n 7^ 1,7V. Expanding the coefficients in plane 
waves a, (3 cx e 4pr+lPzZ we find the energy in terms of in- 
plane p and transverse momentum p z {d is the interlayer 
distance)^ 
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[l-2/ 7iP ]. (11) 



where 7 labels different states for given p, while / 7 . p is 
the distribution function. 

To calculate the supercurrent we use the same ap- 
proach as in Ref. 14 . Consider A = |A|e lkr . Separating 
the order-parameter phase, we put u n — u„(p)e^ p+k / 2 ^ r , 
while v n — v n (p)e i ( p ~ k / 2 ) r . For large N ^> 1 the most 
important corrections come from (p±k/2) Ar . (The exact 
condition for N will be established later.) We have 



K = QivW + (-(p)t- 2 (f 3 E + fi)v F ( Px - ip y )A-, (12) 
C(p)A- + (+(p)t- 2 (f 3 E + fL)v F ( Px + ip y )A+, (13) 
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where p — ^Jp 2 , +p 2 and e l<t> = (p x + ip y )/p. Equations 
Q and Q determine the coefficients^ 1 ^ 



where p = p + f 3 k/2. Equations ((9J, (|T0|) at the outer- 
most layers give 



a+ = Q(p)A + + (-(p)t- 2 (f 3 E + P)v F (p x - iPy)A-,{&) 
&n = (u(p)A~ + &{p)t~ 2 (viE + p)v F ( Px + i Py )A + ,(7) 
where the basis functions are 
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(14) 
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Here 



g = t [v F (p x T iPy)/t] N = e* iN % , C P = t{v FP /t) N 



Cn (P) = [ v f(Px +ip y )/t 
C(P) = [vF(Px-iPy)/t] 
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Using the spinors in the sublattice space, Eqs. (|14j) . 
P^l) can be written as 



Here we include the first-order corrections in energy. 
Having an imaginary momentum p z for v F p < t, these 
solutions decay away from the surfaces and thus they 
describe the surface states. Normalization requires 



Ho + Hi 



where 



tp = Eip , ip 



Ho = f 3 e "''^(TxZp-hfi- 
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[(&+?&+ + (a-)i a'] = 1 
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(16) 

(17) 
(18) 
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In the zero order in k the coefficients satisfy 

A ^ (0) = £ (0) ^ (0) . 

The equation has four solutions 
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Ai 2 e+ m ^ 2 
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Here 



K = V - A) + I A| 2 , = V + fry + | A 



and 
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Normalization is determined by Eq. 
Vpp 2 /t 2 ), the coherence factors are 



I, \c\ 2 =d-\i 



u± 
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U± = 
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The different solutions are orthogonal, 

(^4) = Tr[(^) f 4] = \C\ 2 S lk 



since A^Ai = A\A 3 = 0. The trace is taken over pseudo- 
spin and Nambu indexes. 

If the coefficients A ± are taken in the zero order ap- 
proximation in k, the product a^ct~ in Eq. (|12[) contains 
the exponents e~ 1 ^ and (k x + ik y )e^ 2% ^ and vanishes af- 
ter integration over the momentum directions. Therefore, 
the basis functions £n can be taken in zero approxima- 
tion in k but the coefficients A^ need to be calculated 
up to the first order terms in k. 

The corrections due to the condensate momentum can 
be written as 



(22) 



Equations (JTHJ) - © give 



SE n = , B a p 



\C\ 2 

Corrections to energies are 



\C\ 2 (E a - E P ) 



(23) 
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which is the usual normal-state Doppler shift. 
B12 = B 2 i = B 34 = B 43 = 0, while B 13 = B 31 - 
—B 42 and B 23 = B 32 = B 44 — B 4 \ where 

i([p x k]z) d£ p (u + v_ + v+u_) 



We have 

= — B 24 = 
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B 23 — 



2p dp (E+-Eo) 
i([p x k]z) d^ p (u+w_ — v+vJ) 



2p dp (E+ + Eq) 

The currents Eq. (|11[) at layer n contains the product 
of CC* = (£ P /wFP)e i(Ar ~ 1)0 which is independent of the 
layer number, i.e., of the distance from the surface, and 
the products [(E±fr)/t]C£*C£ « (vfpA) 2( ™ _1) and i( E ± 
fr)/t]Cn*Cn ^ {vFP/t) 2lyN ~ n ^ which decay as functions 
of the distance from the surfaces. All these terms are 
of the order of E/t. We shall see, however, that it is 
the constant term that gives the main contribution to 
the total current through the sample, I = dJ2^ =1 j(n). 
Using J2 P — J(dp/d£, p )pd{;p/2irh we find for the current 
per unit sample width 

£,p d£,p 1 



I = edNk 
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(24) 



To obtain this expression we had to regularize Eq. 
(jlll) which diverges for large £ p . The regularization is 
described in detail in Ref J^. In brief, we subtract the 
normal current which is obtained from the current oper- 
ator taken at energies much higher than A and T. For 
£ p > A, T one has E£ = £ p — /i, Eq = £ p + /x, u = 1, 
and v = 0. Therefore, the diverging part of Eq. (TTTT) is 

£,p d£,p 1 



= -deNk 



1 



2ttS, 



■|C| a f 



(25) 



This contributes to the normal current which, of course, 
turns to zero in the end. Indeed, for A = when the 
particle and hole channels separate, the corrections to 
A^ simply correspond to the full shift of the momentum 
p — > p ± k/2 in the particle (hole) wave functions. As a 
result, the normal current vanishes after the momentum 
integration over the entire Brillouin zonei^. After sub- 
tracting the zero normal current, we arrive at Eq. (|2"4"|) . 
For low temperature T |A|, the last two lines in Eq. 
turn to zero. The total current thus becomes 

"1 2(u + it_ — v + V-) 2 ~ 



deNk. 



£p d£,p 
2nh 
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(26) 



which is similar to the result obtained in RefJ^. 
fi = we have 



For 



deA^k 



2-Kh 
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For large N one can consider p as a slow function as 
compared to £ p . This is equivalent to the assumption 
that d (1 - v 2 F p 2 /t 2 )] /dp = (1 - v 2 F p 2 /t 2 ) (d( p /d P ) 
i.e., that (l - v 2 F p 2 /t 2 ) > (2/N)(v 2 F p 2 /t 2 ). Since ^ ~ A 
we have 

1 - v 2 F p 2 /t 2 = 1 - (A/t)# = (2/N) ln(t/A) 

which holds for A ^> ln(f/A). Therefore, the above con- 
dition is satisfied within the logarithmic approximation. 
Note that neglecting the terms C« *Cn an d Cn*Cn m ^q. 
(fTT|) that decay away from the surfaces is also legitimate 
within the same logarithmic approximation ln(t/A) 3> 1. 
Integrating by parts and using that the integral is deter- 
mined by £ p ~ A we find 



=JVA 2 k 
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The total current does not depend on the sample thick- 
ness Nd as it should be for the surface superconductiv- 
ity. The critical current is determined by max(k x ) ~ ^ 1 
where the coherence length isH £ ~ h/p-ps = frvp/t, 

I c ~ eA ln(t/A)p FB . 

For nonzero fi we find in the same way as in Ref. 14 

em(t/A)k 
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(27) 



Recall that Eq. ([27]) holds for T < |A|. As distinct 
from the case of intrinsic superconductivity in graphene 
considered in Refs i 14 ' 18 ' 19 , the surface superconductivity 
gap |A| is suppressed by doping 1 ^, such that both |A| 
and T c vanish as \x reaches the critical level p c — 2T c q. 

To conclude, we have calculated the zero-temperature 
supercurrent for the surface superconductivity of a flat- 
band multilayered rhombohedral graphene. The super- 
current is finite despite the absence of dispersion of the 
excitation spectrum. The critical current is proportional 
to the zero-temperature gap, i.e., to the superconducting 
critical temperature and to the size of the flat band in the 
momentum space. Nonzero surface supercurrent can be 
responsible for the small Meissner effect and for the sharp 
drop in resistance seen in experiments on graphit e 20 ' 21 . 
The enhanced superconducting density has been reported 
on twin boundaries in Ba^ei-^Co^Asz^. This ob- 
servation can also be considered as indications towards 
surface superconductivity described by our theory. 
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